;; Define SICP package
(defpackage :sicp
  (:use :common-lisp)
  (:export :sicp))

(in-package :sicp)


;; Exercise 1.1.  Below is a sequence of expressions. What is the result printed by
;; the interpreter in response to each expression? Assume that the sequence is to
;; be evaluated in the order in which it is presented.
10

(+ 5 3 4)

(- 9 1)

(/ 6 2)

(setf *a* 3)

(setf *b* (+ *a* 1))

(+ *a* *b* (* *a* *b*))

(eq *a* *b*)

(if (and (> *b* *a*) (< *b* (* *a* *b*)))
    *b*
    *a*)

(cond ((eq *a* 4) 6)
      ((eq *b* 4) (+ 6 7 *a*))
      (t 25))

(+ 2 (if (> *b* *a*)
         *b* *a*))

(* (cond ((> *a* *b*) *a*)
         ((> *b* *a*) *b*)
         (t -1))
   (+ *a* 1))

;; Exercise 1.2.  Translate the following expression into prefix form
(/ (+ 5 4 (- 2 (- 3 (+ 6 4/5))))
   (* 3 (- 6 2) (- 2 7)))

;; Exercise 1.3.  Define a procedure that takes three numbers as arguments and
;; returns the sum of the squares of the two larger numbers.
(defun square (x)
  (* x x))

(defun max1 (x y)
  (if (> x y)
      x
      y))

(defun min1 (x y)
  (if (< x y)
      x
      y))

(defun sum-of-square (x y)
  (+ (square x) (square y)))

(defun sum-of-two-larger-square (x y z)
  (cond ((< x (min1 y z)) (sum-of-square y z))
        (t (sum-of-square x (max1 y z)))))

(sum-of-two-larger-square 1 2 3)
(sum-of-two-larger-square 1 3 2)
(sum-of-two-larger-square 2 1 3)
(sum-of-two-larger-square 2 3 1)
(sum-of-two-larger-square 3 1 2)
(sum-of-two-larger-square 3 2 1)

;; Exercise 1.4 Exercise 1.4.  Observe that our model of evaluation allows for
;; combinations whose operators are compound expressions. Use this observation
;; to describe the behavior of the following procedure:

;; (define (a-plus-abs-b a b)
;;   ((if (> b 0) + -) a b))
(defun a-plus-abs-b (a b)
  (funcall (if (> b 0) #'+ #'-) a b))

(a-plus-abs-b 1 2)
(a-plus-abs-b 1 -2)


;; Exercise 1.7.  The good-enough? test used in computing square roots will not
;; be very effective for finding the square roots of very small numbers. Also,
;; in real computers, arithmetic operations are almost always performed with
;; limited precision. This makes our test inadequate for very large
;; numbers. Explain these statements, with examples showing how the test fails
;; for small and large numbers. An alternative strategy for implementing
;; good-enough? is to watch how guess changes from one iteration to the next and
;; to stop when the change is a very small fraction of the guess. Design a
;; square-root procedure that uses this kind of end test. Does this work better
;; for small and large numbers?
(defun sqrt-iter (guess x)
  (format t "~f " guess)
  (if (good-enoughp guess x)
      guess
      (sqrt-iter (improve guess x)
                 x)))

(defun improve (guess x)
  (average guess (/ x guess)))

(defun average (x y)
  (/ (+ x y) 2))

(defun good-enoughp (guess x)
  (< (abs (- (square guess) x)) 0.001L0))

(defun sqrt1 (x)
  (sqrt-iter 1.0 x))

(defun good-enough2p (diff)
  (< (abs diff) 0.001L0))


(defun sqrt-iter2 (guess x)
  (format t "~f " guess)
  (if (good-enough2p (- (improve guess x) guess))
      guess
      (sqrt-iter2 (improve guess x)
                  x)))

(defun sqrt2 (x)
  (sqrt-iter2 1.0 x))

;; Exercise 1.9.  Each of the following two procedures defines a method for
;; adding two positive integers in terms of the procedures inc, which increments
;; its argument by 1, and dec, which decrements its argument by 1.

;; (define (+ a b)
;;   (if (= a 0)
;;       b
;;       (inc (+ (dec a) b))))

;; (define (+ a b)
;;   (if (= a 0)
;;       b
;;       (+ (dec a) (inc b))))

;; Using the substitution model, illustrate the process generated by each
;; procedure in evaluating (+ 4 5). Are these processes iterative or recursive?
(defun add1 (a b)
  (if (eq a 0)
      b
      (1+ (add1 (1- a) b))))

(add1 4 5)
(1+ (add1 3 5))
(1+ (1+ (add1 2 5)))
(1+ (1+ (1+ (add1 1 5))))
(1+ (1+ (1+ (1+ (add1 0 5)))))
(1+ (1+ (1+ (1+ 5))))
(1+ (1+ (1+ 6)))
(1+ (1+ 7))
(1+ 8)
9

(defun add2 (a b)
  (if (eq a 0)
      b
      (add2 (1- a) (1+ b))))
(add2 4 5)
(add2 3 6)
(add2 2 7)
(add2 1 8)
(add2 0 9)
9

;; Exercise 1.10.  The following procedure computes a mathematical function
;; called Ackermann's function.

;; (define (A x y)
;;   (cond ((= y 0) 0)
;;         ((= x 0) (* 2 y))
;;         ((= y 1) 2)
;;         (else (A (- x 1)
;;                  (A x (- y 1))))))

;; What are the values of the following expressions?

;; (A 1 10)

;; (A 2 4)

;; (A 3 3)

;; Consider the following procedures, where A is the procedure defined above:

;; (define (f n) (A 0 n))

;; (define (g n) (A 1 n))

;; (define (h n) (A 2 n))

;; (define (k n) (* 5 n n))

;; Give concise mathematical definitions for the functions computed by the
;; procedures f, g, and h for positive integer values of n. For example, (k n)
;; computes 5n2.

(defun A (x y)
  (cond ((eq y 0) 0)
        ((eq x 0) (* 2 y))
        ((eq y 1) 2)
        (t (A (1- x)
              (A x (1- y))))))

;; (A 1 10)
;; (A 0 (A 1 9))
;; (A 0 (A 0 (A 1 8)))
;; (A 0 (A 0 (A 0 (A 1 7))))
;; (A 0 (A 0 (A 0 (A 0 (A 1 6)))))
;; (A 0 (A 0 (A 0 (A 0 (A 0 (A 1 5))))))
;; (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 1 4)))))))
;; (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 1 3))))))))
;; (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 1 2)))))))))
;; (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 1 1))))))))))
;; (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 2)))))))))
;; (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 4))))))))
;; (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 8)))))))
;; (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 16))))))
;; (A 0 (A 0 (A 0 (A 0 (A 0 32)))))
;; (A 0 (A 0 (A 0 (A 0 64))))
;; (A 0 (A 0 (A 0 128)))
;; (A 0 (A 0 256))
;; (A 0 512)
;; 1024

;; (A 2 4)
;; (A 1 (A 2 3))
;; (A 1 (A 1 (A 2 2)))
;; (A 1 (A 1 (A 1 (A 2 1))))
;; (A 1 (A 1 (A 1 2)))
;; (A 1 (A 1 (A 0 (A 1 1))))
;; (A 1 (A 1 (A 0 2)))
;; (A 1 (A 1 4))
;; (A 1 (A 0 (A 1 3)))
;; (A 1 (A 0 (A 0 (A 1 2))))
;; (A 1 (A 0 (A 0 (A 0 (A 1 1)))))
;; (A 1 (A 0 (A 0 (A 0 2))))
;; (A 1 (A 0 (A 0 4)))
;; (A 1 (A 0 8))
;; (A 1 16)
;; (A 0 (A 1 15))
;; (A 0 (A 0 (A 1 14)))
;; (A 0 (A 0 (A 0 (A 1 13))))
;; (A 0 (A 0 (A 0 (A 0 (A 1 12)))))
;; (A 0 (A 0 (A 0 (A 0 (A 0 (A 1 11))))))
;; (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 (A 1 10)))))))
;; (A 0 (A 0 (A 0 (A 0 (A 0 (A 0 1024))))))
;; (A 0 (A 0 (A 0 (A 0 (A 0 2048)))))
;; (A 0 (A 0 (A 0 (A 0 4096))))
;; (A 0 (A 0 (A 0 8192)))
;; (A 0 (A 0 16384))
;; (A 0 32768)
;; 65536

;; (A 3 3)
;; (A 2 (A 3 2))
;; (A 2 (A 2 (A 3 1)))
;; (A 2 (A 2 2))
;; (A 2 (A 1 (A 2 1)))
;; (A 2 (A 1 2))
;; (A 2 (A 0 (A 1 1)))
;; (A 2 (A 0 2))
;; (A 2 4)

;; 2n
(defun A-f (n)
  (A 0 n))                              

;; 2^n
(defun A-g (n)
  (A 1 n))

;; 2_1^2_2^2_3..2_n
(defun A-h (n)
  (A 2 n))

;; Exercise 1.11.  A function f is defined by the rule that f(n) = n if n<3 and
;; f(n) = f(n - 1) + 2f(n - 2) + 3f(n - 3) if n> 3. Write a procedure that
;; computes f by means of a recursive process. Write a procedure that computes f
;; by means of an iterative process.
(defun fb3 (n)
  (if (< n 3)
      n
      (+ (fb3 (- n 1))
         (* 2 (fb3 (- n 2)))
         (* 3 (fb3 (- n 3))))))

(defun fb3-2 (n)
  (fb3-iter 0 1 2 n))

(defun fb3-iter (a b c count)
  (if (eq count 0)
       a
       (fb3-iter b c (+ c (* 2 b) (* 3 a)) (1- count))))

;; Exercise 1.12.  The following pattern of numbers is called Pascal's triangle.

;; The numbers at the edge of the triangle are all 1, and each number inside the
;; triangle is the sum of the two numbers above it.35 Write a procedure that
;; computes elements of Pascal's triangle by means of a recursive process.
(defun pascal-triangle (n k)
  (cond ((eq k 1) 1)
        ((eq k n) 1)
        (t (+ (pascal-triangle (1- n) (1- k))
              (pascal-triangle (1- n) k)))))


(defun pascal-triangle-2 (n k)
  (labels ((adj-add (lst)
             (cond ((null (car lst)) nil)
                   ((null (cdr lst)) '(1))
                   (t (cons (+ (car lst) (cadr lst)) (adj-add (cdr lst))))))
           (iter (acc row)
             (if (= row n)
                 acc
                 (iter (cons 1 (adj-add acc)) (1+ row)))))
    (nth (1- k) (iter '(1) 1))))

;; Exercise 1.16.  Design a procedure that evolves an iterative exponentiation
;; process that uses successive squaring and uses a logarithmic number of steps,
;; as does fast-expt. (Hint: Using the observation that (bn/2)2 = (b2)n/2, keep,
;; along with the exponent n and the base b, an additional state variable a, and
;; define the state transformation in such a way that the product a bn is
;; unchanged from state to state. At the beginning of the process a is taken to
;; be 1, and the answer is given by the value of a at the end of the process. In
;; general, the technique of defining an invariant quantity that remains
;; unchanged from state to state is a powerful way to think about the design of
;; iterative algorithms.)
(defun expr-fast (b n)
  (expr-fast-iter 1 b n))

(defun expr-fast-iter (a b n)
  (if (eq n 0)
      a
      (let ((acc (if (evenp n)
                     a
                     (* a b)))
            (h (if (evenp n)
                   (/ n 2)
                   (/ (1- n) 2))))
        (expr-fast-iter acc (square b) h))))

;; Exercise 1.17.  The exponentiation algorithms in this section are based on
;; performing exponentiation by means of repeated multiplication. In a similar
;; way, one can perform integer multiplication by means of repeated
;; addition. The following multiplication procedure (in which it is assumed that
;; our language can only add, not multiply) is analogous to the expt procedure:

;; (define (* a b)
;;   (if (= b 0)
;;       0
;;       (+ a (* a (- b 1)))))

;; This algorithm takes a number of steps that is linear in b. Now suppose we
;; include, together with addition, operations double, which doubles an integer,
;; and halve, which divides an (even) integer by 2. Using these, design a
;; multiplication procedure analogous to fast-expt that uses a logarithmic
;; number of steps.
(defun multi (a b)
  (if (eq 0 b)
      0
      (+ a (multi a (1- b)))))

(defun double (x)
  (+ x x))

(defun halve (x)
  (/ x 2))

(defun fast-multi (a b)
    (cond ((eq 0 b) 0)
          ((evenp b) (double (fast-multi a (halve b))))
          (t (+ a (fast-multi a (1- b))))))

;; Exercise 1.18.  Using the results of exercises 1.16 and 1.17, devise a
;; procedure that generates an iterative process for multiplying two integers in
;; terms of adding, doubling, and halving and uses a logarithmic number of
;; steps.40
(defun fast-multi-iter (a b)
  (labels ((multi-iter (acc x y)
             (cond ((eq 0 y) acc)
                   ((evenp y) (multi-iter acc (double x) (halve y)))
                   (t (multi-iter (+ acc x) x (1- y))))))
    (funcall #'multi-iter 0 a b)))


;; Exercise 1.19.  There is a clever algorithm for computing the Fibonacci
;; numbers in a logarithmic number of steps. Recall the transformation of the
;; state variables a and b in the fib-iter process of section 1.2.2: a a + b and
;; b a. Call this transformation T, and observe that applying T over and over
;; again n times, starting with 1 and 0, produces the pair Fib(n + 1) and
;; Fib(n). In other words, the Fibonacci numbers are produced by applying Tn,
;; the nth power of the transformation T, starting with the pair (1,0). Now
;; consider T to be the special case of p = 0 and q = 1 in a family of
;; transformations Tpq, where Tpq transforms the pair (a,b) according to a bq +
;; aq + ap and b bp + aq. Show that if we apply such a transformation Tpq twice,
;; the effect is the same as using a single transformation Tp'q' of the same
;; form, and compute p' and q' in terms of p and q. This gives us an explicit
;; way to square these transformations, and thus we can compute Tn using
;; successive squaring, as in the fast-expt procedure. Put this all together to
;; complete the following procedure, which runs in a logarithmic number of
;; steps:41

;; (define (fib n)
;;   (fib-iter 1 0 0 1 n))
;; (define (fib-iter a b p q count)
;;   (cond ((= count 0) b)
;;         ((even? count)
;;          (fib-iter a
;;                    b
;;                    <??>      ; compute p'
;;                    <??>      ; compute q'
;;                    (/ count 2)))
;;         (else (fib-iter (+ (* b q) (* a q) (* a p))
;;                         (+ (* b p) (* a q))
;;                         p
;;                         q
;;                         (- count 1)))))
(defun fast-fib (n)
  (fast-fib-iter 1 0 0 1 n))

(defun fast-fib-iter (a b p q count)
  (cond ((= count 0) b)
        ((evenp count)
         (fast-fib-iter a
                        b
                        ;; p' = p^2 + q^2
                        (+ (square p) (square q))
                        ;; q' = q^2 + 2pq 
                        (+ (square q) (* 2 p q))
                        (/ count 2)))
        (t (fast-fib-iter (+ (* b q) (* a q) (* a p))
                          (+ (* b p) (* a q))
                          p
                          q
                          (1- count)))))

;; Exercise 1.20.  The process that a procedure generates is of course dependent
;; on the rules used by the interpreter. As an example, consider the iterative
;; gcd procedure given above. Suppose we were to interpret this procedure using
;; normal-order evaluation, as discussed in section 1.1.5. (The
;; normal-order-evaluation rule for if is described in exercise 1.5.) Using the
;; substitution method (for normal order), illustrate the process generated in
;; evaluating (gcd 206 40) and indicate the remainder operations that are
;; actually performed. How many remainder operations are actually performed in
;; the normal-order evaluation of (gcd 206 40)? In the applicative-order
;; evaluation?
(defun gcd1 (a b)
  (let ((x (max (abs a) (abs b)))
        (y (min (abs a) (abs b))))
    (if (eq 0 y)
        x
        (gcd1 y (mod x y)))))

;; (gcd1 206 40)
;; ;; Normal-order
;; (gcd1 40 (mod 206 40))                  ; 0 + 1
;; (gcd1 (mod 206 40) (mod 40 (mod 206 40))) ; 1 + 1 = 2
;; (gcd1 (mod 40 (mod 206 40)) (mod (mod 206 40) (mod 40 (mod 206 40)))) ; 2 + 1 + 1 = 4
;; (gcd1 (mod (mod 206 40) (mod 40 (mod 206 40))) (mod (mod 40 (mod 206 40)) (mod (mod 206 40) (mod 40 (mod 206 40))))) ; 4 + 2 + 1 = 7
;; (mod (mod 206 40) (mod 40 (mod 206 40))) ; 4
;; 2
;; ;; (+ 1 2 4 7 4) => 18 times

;; Applicative-order
;; (gcd1 40 (mod 206 40))
;; (gcd1 40 6)
;; (gcd1 6 (mod 40 6))
;; (gcd1 6 4)
;; (gcd1 4 (mod 6 4))
;; (gcd1 4 2)
;; (gcd1 2 (mod 4 2))
;; (gcd1 2 0)
;; 2

;; 4 times.

;; Exercise 1.21.  Use the smallest-divisor procedure to find the smallest
;; divisor of each of the following numbers: 199, 1999, 19999.
(defun smallest-divisor (n)
  (find-divisor n 2))

(defun find-divisor (n test-divisor)
  (cond ((devidesp n test-divisor) test-divisor)
        ((> (square test-divisor) n) n)
        (t (find-divisor n (1+ test-divisor)))))

(defun devidesp (n d)
  (eq 0 (mod n d)))

;; (smallest-divisor 199)
;; (smallest-divisor 1999)
;; (smallest-divisor 19999)

;; Exercise 1.22.  Most Lisp implementations include a primitive called runtime
;; that returns an integer that specifies the amount of time the system has been
;; running (measured, for example, in microseconds). The following
;; timed-prime-test procedure, when called with an integer n, prints n and
;; checks to see if n is prime. If n is prime, the procedure prints three
;; asterisks followed by the amount of time used in performing the test.

;; (define (timed-prime-test n)
;;   (newline)
;;   (display n)
;;   (start-prime-test n (runtime)))
;; (define (start-prime-test n start-time)
;;   (if (prime? n)
;;       (report-prime (- (runtime) start-time))))
;; (define (report-prime elapsed-time)
;;   (display " *** ")
;;   (display elapsed-time))

;; Using this procedure, write a procedure search-for-primes that checks the
;; primality of consecutive odd integers in a specified range. Use your
;; procedure to find the three smallest primes larger than 1000; larger than
;; 10,000; larger than 100,000; larger than 1,000,000. Note the time needed to
;; test each prime. Since the testing algorithm has order of growth of (n), you
;; should expect that testing for primes around 10,000 should take about 10
;; times as long as testing for primes around 1000. Do your timing data bear
;; this out? How well do the data for 100,000 and 1,000,000 support the n
;; prediction? Is your result compatible with the notion that programs on your
;; machine run in time proportional to the number of steps required for the
;; computation?
(defun primep (n)
  (eq (smallest-divisor n) n))

(defun timed-prime-test (p n)
  (format t "~d" n)
  (start-prime-test p n (get-internal-run-time)))

(defun start-prime-test (p n start-time)
  (if (funcall p n)
      (report-prime (- (get-internal-run-time) start-time))))

(defun report-prime (elapsed-time)
  (format t " *** ~d" elapsed-time))

;; (timed-prime-test #'primep 1009)
;; (timed-prime-test #'primep 1013)
;; (time (timed-prime-test #'primep 1019))
;; (timed-prime-test #'primep 10007)
;; (timed-prime-test #'primep 10009)
;; (time (timed-prime-test #'primep 10037))
;; (timed-prime-test #'primep 100003)
;; (timed-prime-test #'primep 100019)
;; (time (timed-prime-test #'primep 100043))
;; (timed-prime-test #'primep 1000003)
;; (timed-prime-test #'primep 1000033)
;; (time (timed-prime-test #'primep 1000037))
;; (timed-prime-test #'primep 10000019)
;; (timed-prime-test #'primep 10000079)
;; (time (timed-prime-test #'primep 10000103))

;; Exercise 1.23.  The smallest-divisor procedure shown at the start of this
;; section does lots of needless testing: After it checks to see if the number
;; is divisible by 2 there is no point in checking to see if it is divisible by
;; any larger even numbers. This suggests that the values used for test-divisor
;; should not be 2, 3, 4, 5, 6, ..., but rather 2, 3, 5, 7, 9, .... To implement
;; this change, define a procedure next that returns 3 if its input is equal to
;; 2 and otherwise returns its input plus 2. Modify the smallest-divisor
;; procedure to use (next test-divisor) instead of (+ test-divisor 1). With
;; timed-prime-test incorporating this modified version of smallest-divisor, run
;; the test for each of the 12 primes found in exercise 1.22. Since this
;; modification halves the number of test steps, you should expect it to run
;; about twice as fast. Is this expectation confirmed? If not, what is the
;; observed ratio of the speeds of the two algorithms, and how do you explain
;; the fact that it is different from 2?
(defun smallest-divisor-fast (n)
  (find-divisor-fast n 2))

(defun next-test-divisor (n)
  (if (eq n 2)
      3
      (+ 2 n)))

(defun find-divisor-fast (n test-divisor)
  (cond ((devidesp n test-divisor) test-divisor)
        ((> (square test-divisor) n) n)
        (t (find-divisor n (next-test-divisor test-divisor)))))

(defun fast-primep (n)
  (eq (smallest-divisor-fast n) n))

;; (timed-prime-test #'fast-primep 1009)
;; (timed-prime-test #'fast-primep 1013)
;; (time (timed-prime-test #'fast-primep 1019))
;; (timed-prime-test #'fast-primep 10007)
;; (timed-prime-test #'fast-primep 10009)
;; (time (timed-prime-test #'fast-primep 10037))
;; (timed-prime-test #'fast-primep 100003)
;; (timed-prime-test #'fast-primep 100019)
;; (time (timed-prime-test #'fast-primep 100043))
;; (timed-prime-test #'fast-primep 1000003)
;; (timed-prime-test #'fast-primep 1000033)
;; (time (timed-prime-test #'fast-primep 1000037))
;; (timed-prime-test #'fast-primep 10000019)
;; (timed-prime-test #'fast-primep 10000079)
;; (time (timed-prime-test #'fast-primep 10000103))

;; Exercise 1.24.  Modify the timed-prime-test procedure of exercise 1.22 to use
;; fast-prime? (the Fermat method), and test each of the 12 primes you found in
;; that exercise. Since the Fermat test has (log n) growth, how would you expect
;; the time to test primes near 1,000,000 to compare with the time needed to
;; test primes near 1000? Do your data bear this out? Can you explain any
;; discrepancy you find?
(defun expmod (base exp m)
  (cond ((eq 0 exp) 1)
        ((evenp exp)
         (mod (square (expmod base (/ exp 2) m)) m))
        (t (mod (* base (expmod base (1- exp) m)) m))))


(defun fermat-test (n)
  (labels ((try-it (x)
             (eq (expmod x n n) x)))
    (funcall #'try-it (random n))))

(defun fermat-primep (n &optional (times 10))
  (cond ((eq 0 times) t)
        ((fermat-test n) (fermat-primep n (1- times)))
        (t nil)))

;; (timed-prime-test #'fermat-primep 1009)
;; (timed-prime-test #'fermat-primep 1013)
;; (time (timed-prime-test #'fermat-primep 1019))
;; (timed-prime-test #'fermat-primep 10007)
;; (timed-prime-test #'fermat-primep 10009)
;; (time (timed-prime-test #'fermat-primep 10037))
;; (timed-prime-test #'fermat-primep 100003)
;; (timed-prime-test #'fermat-primep 100019)
;; (time (timed-prime-test #'fermat-primep 100043))
;; (timed-prime-test #'fermat-primep 1000003)
;; (timed-prime-test #'fermat-primep 1000033)
;; (time (timed-prime-test #'fermat-primep 1000037))
;; (timed-prime-test #'fermat-primep 10000019)
;; (timed-prime-test #'fermat-primep 10000079)
;; (time (timed-prime-test #'fermat-primep 10000103))

;; Exercise 1.25.  Alyssa P. Hacker complains that we went to a lot of extra
;; work in writing expmod. After all, she says, since we already know how to
;; compute exponentials, we could have simply written

;; (define (expmod base exp m)
;;   (remainder (fast-expt base exp) m))

;; Is she correct? Would this procedure serve as well for our fast prime tester?
;; Explain.

;; No, the (fast-expt base exp) is too large to process by computer.

;; Exercise 1.26.  Louis Reasoner is having great difficulty doing exercise
;; 1.24. His fast-prime? test seems to run more slowly than his prime?
;; test. Louis calls his friend Eva Lu Ator over to help. When they examine
;; Louis's code, they find that he has rewritten the expmod procedure to use an
;; explicit multiplication, rather than calling square:

;; (define (expmod base exp m)
;;   (cond ((= exp 0) 1)
;;         ((even? exp)
;;          (remainder (* (expmod base (/ exp 2) m)
;;                        (expmod base (/ exp 2) m))
;;                     m))
;;         (else
;;          (remainder (* base (expmod base (- exp 1) m))
;;                     m))))

;; ``I don't see what difference that could make,'' says Louis. ``I do.'' says
;; Eva. ``By writing the procedure like that, you have transformed the (log n)
;; process into a (n) process.'' Explain.

;; T(n) = log(n) + n = O(n)

;; Exercise 1.27.  Demonstrate that the Carmichael numbers listed in footnote 47
;; really do fool the Fermat test. That is, write a procedure that takes an
;; integer n and tests whether an is congruent to a modulo n for every a<n, and
;; try your procedure on the given Carmichael numbers.
(defun test-carmichael-number (n)
  (labels ((test-iter (a)
             (cond ((eq a n) t)
                   ((eq (expmod a n n) a) (test-iter (1+ a)))
                   (t nil))))
    (test-iter 2)))

;; (test-carmichael-number 561)
;; (test-carmichael-number 1105)
;; (test-carmichael-number 1729)
;; (test-carmichael-number 2465)
;; (test-carmichael-number 2821)
;; (test-carmichael-number 6601)

;; (smallest-divisor 561)
;; (smallest-divisor 1105)
;; (smallest-divisor 1729)
;; (smallest-divisor 2465)
;; (smallest-divisor 2821)
;; (smallest-divisor 6601)

;; Exercise 1.28.  One variant of the Fermat test that cannot be fooled is
;; called the Miller-Rabin test (Miller 1976; Rabin 1980). This starts from an
;; alternate form of Fermat's Little Theorem, which states that if n is a prime
;; number and a is any positive integer less than n, then a raised to the (n -
;; 1)st power is congruent to 1 modulo n. To test the primality of a number n by
;; the Miller-Rabin test, we pick a random number a<n and raise a to the (n -
;; 1)st power modulo n using the expmod procedure. However, whenever we perform
;; the squaring step in expmod, we check to see if we have discovered a
;; ``nontrivial square root of 1 modulo n,'' that is, a number not equal to 1 or
;; n - 1 whose square is equal to 1 modulo n. It is possible to prove that if
;; such a nontrivial square root of 1 exists, then n is not prime. It is also
;; possible to prove that if n is an odd number that is not prime, then, for at
;; least half the numbers a<n, computing an-1 in this way will reveal a
;; nontrivial square root of 1 modulo n. (This is why the Miller-Rabin test
;; cannot be fooled.) Modify the expmod procedure to signal if it discovers a
;; nontrivial square root of 1, and use this to implement the Miller-Rabin test
;; with a procedure analogous to fermat-test. Check your procedure by testing
;; various known primes and non-primes. Hint: One convenient way to make expmod
;; signal is to have it return 0.


(defun nontrivial-square-mod (n p)
  (labels ((test-nontrivial-square (x y)
             (cond ((eq x 1) (mod y p))
                   ((eq x (1- p)) (mod y p))
                   ((eq (mod y p) 1) 0)
                   (t (mod y p)))))
    (test-nontrivial-square n (square n))))

(defun expmod1 (base exp m)
  (cond ((eq 0 exp) 1)
        ((evenp exp)
         (nontrivial-square-mod (expmod1 base (/ exp 2) m) m))
        (t (mod (* base (expmod1 base (1- exp) m)) m))))

(defun miller-rabin-test (n)
  (labels ((try-it (x)
             (eq (expmod1 x (1- n) n) 1)))
    (try-it (random n))))

(defun miller-rabin-primep (n &optional (times 10))
  (cond ((eq 0 times) t)
        ((miller-rabin-test n) (miller-rabin-primep n (1- times)))
        (t nil)))

;; (timed-prime-test #'miller-rabin-primep 561)
;; (timed-prime-test #'miller-rabin-primep 1105)
;; (timed-prime-test #'miller-rabin-primep 1729)
;; (timed-prime-test #'miller-rabin-primep 2465)
;; (timed-prime-test #'miller-rabin-primep 2821)
;; (timed-prime-test #'miller-rabin-primep 6601)
;; 
;; (timed-prime-test #'miller-rabin-primep 1009)
;; (timed-prime-test #'miller-rabin-primep 1013)
;; (time (timed-prime-test #'miller-rabin-primep 1019))
;; (timed-prime-test #'miller-rabin-primep 10007)
;; (timed-prime-test #'miller-rabin-primep 10009)
;; (time (timed-prime-test #'miller-rabin-primep 10037))
;; (timed-prime-test #'miller-rabin-primep 100003)
;; (timed-prime-test #'miller-rabin-primep 100019)
;; (time (timed-prime-test #'miller-rabin-primep 100043))
;; (timed-prime-test #'miller-rabin-primep 1000003)
;; (timed-prime-test #'miller-rabin-primep 1000033)
;; (time (timed-prime-test #'miller-rabin-primep 1000037))
;; (timed-prime-test #'miller-rabin-primep 10000019)
;; (timed-prime-test #'miller-rabin-primep 10000079)
;; (time (timed-prime-test #'miller-rabin-primep 10000103))

;; Exercise 1.29.  Simpson's Rule is a more accurate method of numerical
;; integration than the method illustrated above. Using Simpson's Rule, the
;; integral of a function f between a and b is approximated as

;; where h = (b - a)/n, for some even integer n, and yk = f(a + kh). (Increasing
;; n increases the accuracy of the approximation.) Define a procedure that takes
;; as arguments f, a, b, and n and returns the value of the integral, computed
;; using Simpson's Rule. Use your procedure to integrate cube between 0 and 1
;; (with n = 100 and n = 1000), and compare the results to those of the integral
;; procedure shown above.
(defun cube (x)
  (* x x x))

(defun simpson-integrate (f a b n)
  (labels ((next (k)
             (1+ k))
           (ak (k)
             (cond ((or (eq k 0) (eq k n)) 1)
                   ((evenp k) 2)
                   (t 4)))
           (yk (k h)
             (funcall f (+ a (* k h))))
           (term (fa fy k h)
             (* (funcall fa k) (funcall fy k h)))
           (simpson-sum (h fterm fak fyk fnext)
             (labels ((sum-iter (k)
                        (if (> k n)
                            0
                            (+ (funcall fterm fak fyk k h) (sum-iter (funcall fnext k))))))
               (/ (* h (sum-iter 0)) 3))))
    (simpson-sum (/ (- b a) n) #'term #'ak #'yk #'next)))

;(simpson-integrate #'cube 0 1 1000)

;; Exercise 1.30.  The sum procedure above generates a linear recursion. The
;; procedure can be rewritten so that the sum is performed iteratively. Show how
;; to do this by filling in the missing expressions in the following definition:

;; (define (sum term a next b)
;;   (define (iter a result)
;;     (if (> a b)
;;         result
;;         (iter (next a) (+ (term a) result))))
;;   (iter a 0))

(defun simpson-integrate-tail-recursive (f a b n)
  (labels ((next (k)
             (1+ k))
           (ak (k)
             (cond ((or (eq k 0) (eq k n)) 1)
                   ((evenp k) 2)
                   (t 4)))
           (yk (k h)
             (funcall f (+ a (* k h))))
           (term (fa fy k h)
             (* (funcall fa k) (funcall fy k h)))
           (simpson-sum (h fterm fak fyk fnext)
             (labels ((sum-iter (k acc)
                        (if (> k n)
                            acc
                            (sum-iter (funcall fnext k) (+ acc (funcall fterm fak fyk k h))))))
               (/ (* h (sum-iter 0 0)) 3))))
    (simpson-sum (/ (- b a) n) #'term #'ak #'yk #'next)))

;(simpson-integrate-tail-recursive #'cube 0 1 1000)

;; Exercise 1.31.
;; a.  The sum procedure is only the simplest of a vast number of similar
;; abstractions that can be captured as higher-order procedures.51 Write an
;; analogous procedure called product that returns the product of the values of
;; a function at points over a given range. Show how to define factorial in
;; terms of product. Also use product to compute approximations to using the
;; formula52

;; b.  If your product procedure generates a recursive process, write one that
;; generates an iterative process. If it generates an iterative process, write
;; one that generates a recursive process.

(defun product (fterm a fnext b)
  (labels ((product-iter (n)
             (if (> n b)
                 1
                 (* (funcall fterm n)
                    (product-iter (funcall fnext n))))))
    (product-iter a)))

(defun product-tail-recursive (fterm a fnext b)
  (labels ((product-iter (n acc)
             (if (> n b)
                 acc
                 (product-iter (funcall fnext n)
                               (* (funcall fterm n) acc)))))
    (product-iter a 1)))


(defun factorial (n)
  (labels ((fterm (x)
             x))
    (product #'fterm 1 #'1+ n)))

(defun wallis-product (n)
  (labels ((fterm (k)
            (expt (/ (+ 2 k) (+ 3 k)) (expt -1 k))))
    (product #'fterm 0 #'1+ n)))

;; Exercise 1.32.  a. Show that sum and product (exercise 1.31) are both special
;; cases of a still more general notion called accumulate that combines a
;; collection of terms, using some general accumulation function:

;; (accumulate combiner null-value term a next b)

;; Accumulate takes as arguments the same term and range specifications as sum
;; and product, together with a combiner procedure (of two arguments) that
;; specifies how the current term is to be combined with the accumulation of the
;; preceding terms and a null-value that specifies what base value to use when
;; the terms run out. Write accumulate and show how sum and product can both be
;; defined as simple calls to accumulate.

;; b. If your accumulate procedure generates a recursive process, write one that
;; generates an iterative process. If it generates an iterative process, write
;; one that generates a recursive process.

(defun accumulate (combiner null-value term a next b)
  (labels ((accumulate-iter (n)
             (if (> n b)
                 null-value
                 (funcall combiner (funcall term n)
                          (accumulate-iter (funcall next n))))))
    (accumulate-iter a)))

(defun accumulate-tail-recursive (combiner null-value term a next b)
  (labels ((accumulate-iter (n acc)
             (if (> n b)
                 acc
                 (accumulate-iter (funcall next n)
                                  (funcall combiner (funcall term n) acc)))))
    (accumulate-iter a null-value)))

(defun general-sum (n)
  (accumulate #'+ 0 #'(lambda (x) x) 1 #'1+ n))

(defun general-wallis-product (n)
  (labels ((fterm (k)
            (expt (/ (+ 2 k) (+ 3 k)) (expt -1 k))))
    (accumulate #'* 1 #'fterm 0 #'1+ n)))

;; Exercise 1.33.  You can obtain an even more general version of accumulate
;; (exercise 1.32) by introducing the notion of a filter on the terms to be
;; combined. That is, combine only those terms derived from values in the range
;; that satisfy a specified condition. The resulting filtered-accumulate
;; abstraction takes the same arguments as accumulate, together with an
;; additional predicate of one argument that specifies the filter. Write
;; filtered-accumulate as a procedure. Show how to express the following using
;; filtered-accumulate:

;; a. the sum of the squares of the prime numbers in the interval a to b
;; (assuming that you have a prime? predicate already written)

;; b. the product of all the positive integers less than n that are relatively
;; prime to n (i.e., all positive integers i < n such that GCD(i,n) = 1).
(defun filtered-accumulate (combiner filterp null-value term a next b)
  (labels ((accumulate-iter (n)
             (if (> n b)
                 null-value
                 (funcall combiner
                          (if (funcall filterp n)
                              (funcall term n)
                              null-value)
                          (accumulate-iter (funcall next n))))))
    (accumulate-iter a)))

(defun filtered-accumulate-tail-recursive (combiner filterp null-value term a next b)
  (labels ((accumulate-iter (n acc)
             (if (> n b)
                 acc
                 (accumulate-iter (funcall next n)
                                  (funcall combiner
                                           (if (funcall filterp n)
                                               (funcall term n)
                                               null-value) acc)))))
    (accumulate-iter a null-value)))

(defun more-general-sum (n)
  (filtered-accumulate #'+ #'(lambda (x) t) 0 #'(lambda (x) x) 1 #'1+ n))

(defun more-general-wallis-product (n)
  (labels ((fterm (k)
            (expt (/ (+ 2 k) (+ 3 k)) (expt -1 k))))
    (filtered-accumulate-tail-recursive #'* #'(lambda (x) t) 1 #'fterm 0 #'1+ n)))

;(more-general-sum 100)
;(time (* 4.0 (more-general-wallis-product 25000)))
(defun sum-of-primer-square (a b)
  (filtered-accumulate #'+ #'primep 0 #'square a #'1+ b))

(defun products-of-relative-primer-to (n)
  (filtered-accumulate #'* #'(lambda (x)
                               (eq (gcd1 n x) 1))
                       1 #'(lambda (x) x) 1 #'1+ (1- n)))

; (eq (factorial (1- 17)) (products-of-relative-primer-to 17))

;; Exercise 1.34.  Suppose we define the procedure

;; (define (f g)
;;   (g 2))

;; Then we have

;; (f square)
;; 4

;; (f (lambda (z) (* z (+ z 1))))
;; 6

;; What happens if we (perversely) ask the interpreter to evaluate the
;; combination (f f)? Explain.
(defun f-1-34 (g)
  (funcall g 2))

; (f-1-34 #'square)
; (f-1-34 #'(lambda (z) (* z (+ z 1))))
; (f-1-34 #'f-1-34) => (funcall #'f-1-34 2) => (funcall 2 2)

;; *Exercise 1.35:* Show that the golden ratio [phi] (section *note 1-2-2::) is
;; *a fixed point of the transformation x |-> 1 + 1/x, and use this fact to
;; *compute [phi] by means of the `fixed-point' procedure.
(defparameter *tolerance* 0.00001d0)

(defun fixed-point (f first-guess)
  (labels ((close-enoughp (v1 v2)
             (< (abs (- v1 v2)) *tolerance*))
           (try (guess)
             (let ((next (funcall f guess)))
               (if (close-enoughp guess next)
                   next
                   (try next)))))
    (try first-guess)))

;; (fixed-point #'(lambda (x) (+ 1 (/ 1 x))) 2.71828d0)

;; Exercise 1.36: Modify `fixed-point' so that it prints the sequence of
;; approximations it generates, using the `newline' and `display' primitives
;; shown in *note Exercise 1-22::.  Then find a solution to x^x = 1000 by
;; finding a fixed point of x |-> `log'(1000)/`log'(x).  (Use Scheme's primitive
;; `log' procedure, which computes natural logarithms.)  Compare the number of
;; steps this takes with and without average damping.  (Note that you cannot
;; start `fixed-point' with a guess of 1, as this would cause division by
;; `log'(1) = 0.)
(defun traced-fixed-point (f first-guess)
  (labels ((close-enoughp (v1 v2)
             (< (abs (- v1 v2)) *tolerance*))
           (try (guess)
             (let ((next (funcall f guess)))
               (format t "~f" guess)
               (fresh-line)
               (if (close-enoughp guess next)
                   next
                   (try next)))))
    (try first-guess)))

; (traced-fixed-point #'(lambda (x) (/ (log 1000) (log x))) pi)
; (traced-fixed-point #'(lambda (x) (average x (/ (log 1000) (log x)))) pi)
; (traced-fixed-point #'(lambda (x) (+ 1 (/ 1 x))) -1.5)

;; *Exercise 1.37:* a. An infinite "continued fraction" is an expression of the
;; form

;;                           N_1
;;                f = ---------------------
;;                               N_2
;;                    D_1 + ---------------
;;                                   N_3
;;                          D_2 + ---------
;;                                D_3 + ...

;; As an example, one can show that the infinite continued fraction expansion
;; with the n_i and the D_i all equal to 1 produces 1/[phi], where [phi] is the
;; golden ratio (described in section *note 1-2-2::).  One way to approximate an
;; infinite continued fraction is to truncate the expansion after a given number
;; of terms.  Such a truncation--a so-called finite continued fraction "k-term
;; finite continued fraction"--has the form

;;             N_1
;;      -----------------
;;                N_2
;;      D_1 + -----------
;;            ...    N_K
;;                + -----
;;                   D_K

;; Suppose that `n' and `d' are procedures of one argument (the term index i)
;; that return the n_i and D_i of the terms of the continued fraction.  Define a
;; procedure `cont-frac' such that evaluating `(cont-frac n d k)' computes the
;; value of the k-term finite continued fraction.  Check your procedure by
;; approximating 1/[phi] using

;;      (cont-frac (lambda (i) 1.0)
;;                 (lambda (i) 1.0)
;;                 k)

;; for successive values of `k'.  How large must you make `k' in order to get an
;; approximation that is accurate to 4 decimal places?

;; b. If your `cont-frac' procedure generates a recursive process, write one
;; that generates an iterative process.  If it generates an iterative process,
;; write one that generates a recursive process.
(defun cont-frac (n d k)
  (labels ((iter (x)
             (let ((term-n (funcall n x))
                   (term-d (funcall d x)))
               (if (eq x k)
                   (/ term-n term-d)
                   (/ term-n (+ term-d (iter (1+ x))))))))
    (iter 1)))

(defun cont-frac-iterate (n d k)
  (labels ((iter (x acc)
             (let ((term-n (funcall n x))
                   (term-d (funcall d x)))
               (if (eq x 1)
                   (/ term-n (+ term-d acc))
                   (iter (1- x) (/ term-n (+ term-d acc)))))))
    (iter k 0)))

;; (expt (/ (+ 1 (sqrt 5.0L0)) 2) -1)
;; (eql (cont-frac #'(lambda (i) 1.0D0) #'(lambda (i) 1.0D0) 17)
;;      (cont-frac-iterate #'(lambda (i) 1.0L0) #'(lambda (i) 1.0L0) 17))

;; *Exercise 1.38:* In 1737, the Swiss mathematician Leonhard Euler published a
;; memoir `De Fractionibus Continuis', which included a continued fraction
;; expansion for e - 2, where e is the base of the natural logarithms.  In this
;; fraction, the n_i are all 1, and the D_i are successively 1, 2, 1, 1, 4, 1,
;; 1, 6, 1, 1, 8, ....  Write a program that uses your `cont-frac' procedure
;; from *note Exercise 1-37:: to approximate e, based on Euler's expansion.

(defun de-fractionibus-continuis (k)
  (+ 2 (cont-frac #'(lambda (i) 1.0d0)
                  #'(lambda (i)
                      (let ((j (- i 2)))
                        (cond ((< i 3) i)
                              ((eq (mod j 3) 0)
                               (+ 2 (* 2 (/ j 3))))
                               (t 1))))
                      k)))

; (de-fractionibus-continuis 10)

;; *Exercise 1.39:* A continued fraction representation of the tangent function
;; was published in 1770 by the German mathematician J.H. Lambert:

;;                    x
;;      tan x = ---------------
;;                      x^2
;;              1 - -----------
;;                        x^2
;;                  3 - -------
;;                      5 - ...

;; where x is in radians.  Define a procedure `(tan-cf x k)' that computes an
;; approximation to the tangent function based on Lambert's formula.  `K'
;; specifies the number of terms to compute, as in *note Exercise 1-37::.
(defun tan-cf (x k)
  (let ((neg-square-x (* -1 (square x))))
    (cont-frac-iterate #'(lambda (i)
                           (if (eq i 1)
                               x
                             neg-square-x))
                       #'(lambda (i)
                           (+ 1.0d0 (* 2 (1- i))))
                       k)))

;; (tan-cf (* pi 7/4) 19)

;; *Exercise 1.40:* Define a procedure `cubic' that can be used together with
;; the `newtons-method' procedure in expressions of the form

;;      (newtons-method (cubic a b c) 1)

;; to approximate zeros of the cubic x^3 + ax^2 + bx + c.
(defparameter *dx* 0.00001d0)

(defun deriv (g)
  (lambda (x)
    (/ (- (funcall g (+ x *dx*)) (funcall g x))
       *dx*)))

;; (funcall (deriv #'cube) 5)
(defun newton-transform (g)
  (lambda (x)
    (- x (/ (funcall g x) (funcall (funcall #'deriv g) x)))))

(defun newtons-method (g guess)
  (fixed-point (newton-transform g) guess))

(defun cubic (a b c)
  (lambda (x) (+ (* a (expt x 3)) (* b (expt x 2)) c)))

;; (funcall (cubic 3 7 11) (newtons-method (cubic 3 7 11) 1))

;; *Exercise 1.41:* Define a procedure `double' that takes a procedure of one
;; argument as argument and returns a procedure that applies the original
;; procedure twice.  For example, if `inc' is a procedure that adds 1 to its
;; argument, then `(double inc)' should be a procedure that adds 2.  What value
;; is returned by

;;      (((double (double double)) inc) 5)

(defun double-proc (g)
  (lambda (x) (funcall g (funcall g x))))

;; (funcall (double-proc #'1+) 0)
;; (funcall (funcall (double-proc #'double-proc) #'1+) 0)
;; (funcall (funcall (double-proc (double-proc #'double-proc)) #'1+) 0)

;; *Exercise 1.42:* Let f and g be two one-argument functions.  The
;; "composition" f after g is defined to be the function x |-> f(g(x)).  Define
;; a procedure `compose' that implements composition.  For example, if `inc' is
;; a procedure that adds 1 to its argument,

;;      ((compose square inc) 6)
;;      49

(defun compose-proc (f g)
  (lambda (x) (funcall f (funcall g x))))

;; (funcall (compose-proc #'square #'1+) 6)

;; *Exercise 1.43:* If f is a numerical function and n is a positive integer,
;; then we can form the nth repeated application of f, which is defined to be
;; the function whose value at x is f(f(...(f(x))...)).  For example, if f is
;; the function x |-> x + 1, then the nth repeated application of f is the
;; function x |-> x + n.  If f is the operation of squaring a number, then the
;; nth repeated application of f is the function that raises its argument to the
;; 2^nth power.  Write a procedure that takes as inputs a procedure that
;; computes f and a positive integer n and returns the procedure that computes
;; the nth repeated application of f.  Your procedure should be able to be used
;; as follows:

;;      ((repeated square 2) 5)
;;      625

;; Hint: You may find it convenient to use `compose' from *note Exercise 1-42::.
(defun repeated-proc (f n)
  (labels ((iter (g i)
             (if (eq i n)
                 g
                 (iter (compose-proc f g) (1+ i)))))
    (iter f 1)))

;; (funcall (repeated-proc #'square 2) 5)

;; *Exercise 1.44:* The idea of "smoothing" a function is an important concept
;; in signal processing.  If f is a function and dx is some small number, then
;; the smoothed version of f is the function whose value at a point x is the
;; average of f(x - dx), f(x), and f(x + dx).  Write a procedure `smooth' that
;; takes as input a procedure that computes f and returns a procedure that
;; computes the smoothed f.  It is sometimes valuable to repeatedly smooth a
;; function (that is, smooth the smoothed function, and so on) to obtained the
;; "n-fold smoothed function".  Show how to generate the n-fold smoothed
;; function of any given function using `smooth' and `repeated' from *note
;; Exercise 1-43::.
(defun smooth (f)
  (lambda (x)
    (/ (+ (funcall f (+ x *dx*)) (funcall f x) (funcall f (- x *dx*))) 3)))

(defun n-fold-smooth (f n)
    (funcall (repeated-proc #'smooth n) f))

;(funcall (n-fold-smooth #'sin 3) 0.5d0)

;; *Exercise 1.45:* We saw in section *note 1-3-3:: that attempting to compute
;; square roots by naively finding a fixed point of y |-> x/y does not converge,
;; and that this can be fixed by average damping.  The same method works for
;; finding cube roots as fixed points of the average-damped y |-> x/y^2.
;; Unfortunately, the process does not work for fourth roots--a single average
;; damp is not enough to make a fixed-point search for y |-> x/y^3 converge.  On
;; the other hand, if we average damp twice (i.e., use the average damp of the
;; average damp of y |-> x/y^3) the fixed-point search does converge.  Do some
;; experiments to determine how many average damps are required to compute nth
;; roots as a fixed-point search based upon repeated average damping of y |->
;; x/y^(n-1).  Use this to implement a simple procedure for computing nth roots
;; using `fixed-point', `average-damp', and the `repeated' procedure of
;; *notExercise 1-43 Assume that any arithmetic operations you need are
;; available as primitives.
(defun average-damp (f)
  (lambda (x)
    (average (funcall f x) x)))

(defun nth-root (x n)
  (labels ((fy (y)
             (/ x  (expt y (1- n)))))
    ;; repeat i times, 2^i >= n
    (fixed-point (funcall (repeated-proc #'average-damp (ceiling (log n 2)))
                          #'fy)
                 1.0d0)))

;; *Exercise 1.46:* Several of the numerical methods described in this chapter
;; are instances of an extremely general computational strategy known as
;; "iterative improvement".  Iterative improvement says that, to compute
;; something, we start with an initial guess for the answer, test if the guess
;; is good enough, and otherwise improve the guess and continue the process
;; using the improved guess as the new guess.  Write a procedure
;; `iterative-improve' that takes two procedures as arguments: a method for
;; telling whether a guess is good enough and a method for improving a guess.
;; `Iterative-improve' should return as its value a procedure that takes a guess
;; as argument and keeps improving the guess until it is good enough.  Rewrite
;; the `sqrt' procedure of section *note 1-1-7:: and the `fixed-point' procedure
;; of section *note 1-3-3:: in terms of `iterative-improve'.
(defun interactive-improve (nicep improve guess)
  (labels ((iter (v)
             (if (funcall nicep v)
                 v
                 (iter (funcall improve v)))))
    (iter guess)))

(defun general-fixed-point (f guess)
  (interactive-improve #'(lambda (v1 v2)
                           (< (abs (- v1 v2)) *tolerance*))
                       f guess))

(defun general-sqrt (n)
  (interactive-improve #'(lambda (x)
                           (< (abs (- (square x) n)) 1e-8))
                       #'(lambda (x) (average x (/ n x)))
                       1.0L0))